Generation of dynamic structures in nonequilibrium reactive bilayers 
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We present a nonequlibrium approach for the study of a flexible bilayer whose two components 
induce distinct curvatures. In turn, the two components are interconverted by an externally pro- 
moted reaction. Phase separation of the two species in the surface results in the growth of domains 
characterized by different local composition and curvature modulations. This domain growth is 
limited by the effective mixing due to the interconversion reaction, leading to a finite characteristic 
domain size. In addition to these effects, first introduced in our earlier work [Phys. Rev. E 71, 
051906 (2005)], the important new feature is the assumption that the reactive process actively affects 
the local curvature of the bilayer. Specifically, we suggest that a force energetically activated by 
external sources causes a modification of the shape of the membrane at the reaction site. Our results 
show the appearance of a rich and robust dynamical phenomenology that includes the generation of 
traveling and/or oscillatory patterns. Linear stability analysis, amplitude equations and numerical 
simulations of the model kinetic equations confirm the occurrence of these spatiotemporal behaviors 
in nonequilibrium reactive bilayers. 



T3 

C 

O 
o 



> 

m 
o 



in 
o 



T3 
C 

o 
o 



PACS numbers: 87.16.Dg, 07.05. Tp, 82.45.Mp 



I. INTRODUCTION 

Due to hydrophobic repulsions, amphiphilic molecules 
such as lipids spontaneously aggregate in water to form 
bilayer membranes. These bilayers typically exhibit an 
in-plane fluidlike nature, and are highly flexible surfaces 
so they can display a large variety of conformations. The 
experimental study of equilibrium lipid bilayers has at- 
tracted a great deal of attention in the past decades 0,0- 
Issues such as membrane conformational behavior, shape 
fluctuations, fusion and fission, and phase segregation in 
multicomponent bilayers have been extensively investi- 
gated [l], |2, |3( ■ Among many other manipulation tech- 
niques, the fabrication of synthetic giant vesicles and 
planar lipid membranes Q as well as micropipet aspira- 
tion p| have been used to perform these studies, leading 
to a fairly detailed current understanding of equilibrium 
membranes. 

Parallel to the experimental advances, theoretical 
modeling of equilibrium flexible membranes has also 
come a long way. Typical theoretical approaches are 
based on the search for minimum energy conformations 
according to the Canham-Helfrich bending elasticity de- 
scription 0, 0] of a tensionless membrane, plus (option- 
ally) the corresponding area-difference contributions for 
closed surfaces 8] . More recently, it has been recognized 
that internal degrees of freedom such as chemical compo- 
sition can crucially influence the shape of the bilayer, es- 
pecially when dealing with multicomponent membranes 
undergoing phase separation. In particular, models con- 
taining a local coupling between composition and curva- 
ture have been developed 0, EIH lH El H, result- 



ing in an interesting outcome, namely, the appearance of 
phase separating domains that grow and adopt distinct 
membrane curvatures. This is understood as the initi- 
ation mechanism for budding phenomena in real mem- 
branes. Kinetic mesoscopic schemes, Monte Carlo ap- 
proaches, and a variety of discrete dynamic algorithms 
have been proposed on this basis, and have mainly been 
devoted to the study of the equilibrium behavior and the 
dynamics toward equilibrium of flexible bilayers. 

The ultimate motivation for the study of lipid bilay- 
ers lies in their likeness to cell membranes. A biolog- 
ical membrane is a complex mixture of lipids, sterols 
and proteins that represents the main structural com- 
ponent of the cell architecture. Rather than an inert 
static boundary, a membrane has to be viewed as a dy- 
namical surface directly and actively involved in many 
biological cell processes .15]. In addition to their compo- 
sitional complexity, cell membranes are also continuously 
subjected to nonequilibrium driving forces, chemical gra- 
dients, and energy fluxes, so it should be recognized that 
nonequilibrium behavior may underlie many aspects of 
their dynamics [Til IT?! ITsf . Therefore, although experi- 
mental and theoretical studies on thermally equilibrated 
membranes have had considerable success and are a good 
starting point, nonequilibrium approaches are fundamen- 
tal for a more accurate understanding of the dynamical 
properties of both natural membranes and synthetic lipid 
bilayers. 

Excellent work on nonequilibrium lipid bilayers has 
recently been presented by the experimental groups of 
J. Prost and H.-G. Dobereiner. J. Prost et al. have 
studied the effects of active proteins inserted in a vesi- 
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cle membrane. Acting as ion pumps externally acti- 
vated by light, these proteins provide a nonthermal en- 
ergy source that directly affects membrane shape be- 
havior 1 1 9| . As an example of another nonequilibrium 
source, H.-G. Dobereiner et al. investigated the mor- 
phological transformations of a vesicle membrane due to 
a photoinduced chemical reaction that modifies the lo- 
cal lipid composition of the bilayer |20| . Some progress 
has also been made on the modeling front. For example, 
parallel to their experimental findings, a novel nonequi- 
librium scheme has been proposed by J. Prost et al. 
that successf ully exp lains the observed experimental phe- 
nomenology |U|22|. In addition to this proposal, other 
nonequilibrium situations have been modeled very re- 
cently, such as those of membranes confined between par- 
allel plates |23|. membranes with active inclusions |24| . 
and bilayers near repulsive walls plf . 

We are mainly interested in the nonequilibrium situa- 
tion proposed by H.-G. Dobereiner et al., where the bi- 
layer's own lipid constituents are chemically transformed, 
leading to changes in the membrane curvature. The 
influence of lipid composition on the curvature of cell 
membranes has recently been established in cell fusion 
processes where high-curvature lipids play an important 
role (2(| . Chemically induced lipid modifications are also 
known to be responsible for membrane shape transforma- 
tions involved in nervous synaptic processes such as the 
formation of microvesicles that release the neurotrans- 
mitter to the gap between two nerve cells [27I. Fur- 
thermore, in addition to biological membranes, an un- 
derstanding of the coupling between chemical reactions 
and the interfacial curvature is essential in elucidating 
dynamical processes in artificial bilayers that may be use- 
ful in nanotechnology applications. As an example, self- 
assembly techniques based on biomineralization can lead 
to hierarchically organized materials built from reactive 
amphiphilic interfaces [28[ . 

Based on the nonequilibrium ingredients described 
above, in a previous paper p9j we presented a model 
for a two-component reactive membrane that exhibits 
a particular compositional and morphological organiza- 
tion. In that model the two lipid constituents have a 
distinctly different shape, so that they are able to pro- 
duce distinct curvatures in the membrane. Moreover, 
they are assumed to be thermodynamically immiscible, 
so they spontaneously induce the development of phase 
separating domains. Additionally, an externally induced 
reaction (promoted by a nonequilibrium source such as 
a chemical flux, an applied light, etc.) interconverts the 
two lipids. The combination of these ingredients leads to 
stationary patterns involving heterogeneous modulations 
of composition and curvature. The novelty of the re- 
sulting pattern organization is that although stationary, 
the patterns are generated in a nonequilibrium context, 
that is, they are actively maintained by the competition 
between thermodynamic (equilibrium) phase separation 
and the environmentally induced (nonequilibrium) reac- 
tion. At thermal equilibrium, the low affinity between 



the two immiscible components generates domains with 
a characteristic composition and curvature. In the ab- 
sence of reaction, these domains grow indefinitely and co- 
alesce until complete segregation into two large domains 
is achieved. However, the reactive process converts one 
lipid component into the other, resulting in a large-scale 
mixing mechanism that halts the growth of the segre- 
gated structures when the mixing action balances the 
short-scale ordering effect of phase separation. This bal- 
ance leads to a stationary state with patterns of a finite 
size. 

In our previous model 29], we assumed that the non- 
thermal energy contribution that promoted the lipid in- 
terconversion reaction simply dissipates to the medium. 
In this paper we generalize our original model to the case 
where such an energetic process directly alters the mem- 
brane shape by exerting local forces on it, that is, the 
induced reactive process locally 'kicks' the membrane. 
Conformational changes of the membrane constituents 
due to the reaction can reasonably be expected to have 
mechanical effects on the membrane. We model this ef- 
fect in a very generic way, and will demonstrate that 
the inclusion of this new ingredient is essential for the 
generation of dynamic spatial organization. As a typi- 
cal feature of soft-matter systems, robust spatiotemporal 
phenomena such as those displayed in this paper are in- 
teresting in their applicability to real bilayer membranes. 

Although our ideas are inspired by phenomena ob- 
served in real cellular systems, we hasten to disclaim a 
direct analogy between our model and in vivo systems 
which are far more complex than our stripped model. 
We only wish to suggest that the mechanisms proposed 
herein might play a role in real membranes in which 
the energy sources for the local membrane-shape-altering 
processes might include ATP consumption, the asborp- 
tion of light, or any of a number of other possible energy 
providing mechanisms. We point out here and again later 
in the paper that the spatiotemporal structures of inter- 
est require parameters that constrain the range of appli- 
cability of the model. In particular, we will see that elas- 
tic properties (bending rigidity) play a central role, and 
that the most likely candidates for the observation of our 
predictions due to their flexibility are in vitro lipid bilay- 
ers formed from single- and/or short-chain lipid surfac- 
tants, especially those composed of surfactants of rather 
distinct shapes that are interconverted by a reactive pro- 
cess that can be controlled experimentally. We discuss 
these and other possibilities for experimental realization 
later in the paper. 

The outline of this paper is as follows. Section [H] 
presents our basic mesoscopic free energy description as 
well as the derivation of the kinetic equations for the com- 
position and curvature variables in the proposed model. 
The linear stability analysis and the amplitude equation 
analysis are presented in Sec. II I II with the details of the 
derivation of the latter relegated to the Appendix. Or- 
ganized spatiotemporal behaviors are predicted in some 
regions of parameter space provided that the parame- 
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ter describing the activity of the reactive process on the 
membrane dynamics is sufficiently large. In good agree- 
ment with the analytical predictions, numerical simula- 
tions exhibiting traveling and/or oscillating domains are 
shown in Sec. IIVI Finally, in Sec. we summarize the 
main conclusions of this paper. 



II. THE MODEL 

In our earlier work we considered the simplest scenario 
where one layer of the membrane (say the outer one) 
was composed of two differently shaped lipids A and B, 
whereas the other layer was composed of a single compo- 
nent without any curvature effect. We also disregarded 
any lipid flip-flop exchange between the inner and outer 
layers. Under these assumptions, we modeled the prop- 
erties of the non-active membrane as a two-dimensional 
surface characterized by the local concentration differ- 
ence <f) between the A and B components and the local 
extrinsic curvature H of the surface. At this point we 
do not specifically identify the components A and B in 
any further detail since this is to be seen as a schematic 
and highly streamlined representation of any number of 
possible realizations. Indeed, A and B need not even 
be single compounds; for example, A may be a group of 
lipids or biomolccules some of which induce a particular 
curvature, while B may be another group that induces a 
different curvature. A and B could be the isomers of an 
amphiphillic azobenzene derivative, or groups of lipids as 
in Fig. 2 in [^3 or those suggested in [2(|, or simply a 
given lipid with its polar head more or less charged as in 
the experiments in 20]. 

The free energy functional in terms of the two order 
parameters was assumed to be 



dxdy 



|0 4 + ||V0| 2 



(1) 



The first three terms correspond to the typical Ginzburg- 
Landau expansion that leads to phase separation when 
a, /?, and 7 are all positive, with an equilibrium concen- 
tration difference <fi eq — ±y // aj]3 and a typical interface 
length £ = y/~f/a. When a is negative the components 
are completely miscible, leading to a homogeneous state. 
The last term is the elastic energy contribution due to 
the rigidity of the membrane, and k is the bending rigid- 
ity modulus. For simplicity we have assumed that the 
spontaneous (equilibrium) curvature H sp (</>), which re- 
flects the shape asymmetry between the two lipid com- 
ponents, is linear in the order parameter, H sp = 4>Hq. 
In the Monge parametrization [3(1 ] a deformable surface 
is described by (x,y,h(x,y)), where h(x,y) is the dis- 
placement (height) field for the local separation from the 
flat conformation. This representation is valid for sur- 
faces that are nearly flat with only gradual variations 
of h, and allows the approximation H « V h that we 



have incorporated in the free energy functional. For self- 
assembled free membranes, the surface tension contribu- 
tion (^|V/i| 2 ) in the free energy can be neglected, and 
we have not included it in Eq. (|TJ. 

The kinetics of <f> and h are assumed to be driven by the 
free energy functional, and by the chemical reaction that 
interconverts the two lipid components, A ^ B. The 
effect of the reaction is two-fold: it affects the concentra- 
tion difference order parameter </>, and it also affects the 
shape of the membrane. This last effect is our new con- 
tribution in this paper. The kinetics is thus a generalized 
version of that introduced in 123 : 



at 



= DV 2 



dh . b"T _,, , , . . 

at = - A flr + r W-*>K 



r(0-0 o ), (2a) 
(2b) 



The concentration difference 4> follows a conserved 
scheme |3l| with diffusion coefficient D 1 augmented by 
the reaction contributions. The rate parameter T = 
k + + fc_ and the stationary concentration difference pa- 
rameter <p = {k_ — (fc + + k_ ) are determined by the 
forward and backward reaction rate constants k + and 
k-. Note that we have assumed local overdamped re- 
laxational dynamics for the membrane height, that is, 
Rouse-like dynamics, where A denotes a mobility param- 
eter proportional to the inverse of the typical relaxatio 
time Th- This approximation is only valid when the mem- 
brane is immersed in a high viscosity medium and/or 
when the membrane is highly permeable |32|. When hy- 
drodynamic effects can not be neglected, a more realistic, 
albeit complex, description of the membrane dynamics 
(Zimm-like dynamics) is required, where inertial effects 
and a space dependent relaxation parameter are in or- 
der m. 

The last term in Eq. (|2b|) is based on the following hy- 
pothesis. The reaction that interconverts A and B can 
be understood as an isomerization-likc chemical transfor- 
mation, involving a strong modification in the shape of 
the membrane constituents. This process implies the dis- 
placement of parts of these molecules that could have a 
mechanical effect on the local membrane shape. If addi- 
tionally the process is strongly energetically activated by 
an external energy source, it might exert a local force on 
the membrane. A simplifying approximation is to assume 
that the forward and backward reactions exert opposite 
forces on the membrane. These forces are assumed to act 
locally for a negligible period of time (the time needed 
to complete the reaction is much shorter than any other 
time scale of the system), in the same direction as the 
preferred curvature of the reaction product component, 
that is, positive (outwards) for A —> B, and negative 
(inwards) for B — > A, see Fig. ^ This is modeled in a 
generic way |34j by the last term in the membrane height 
kinetic equation, where £ is a parameter accounting for 
the strength of the effect of the reaction on the shape of 
the membrane. Here we take the parameter £ to have the 
same sign as Hq (in Fig.^both are positive). Active pro- 
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teins are known to act as force centers when inserted in 
lipid bilayers 0], an d other experimental studies [Tillill 
also provide evidence that reactive processes may locally 
modify the membrane shape in red blood cells. How- 
ever, we must recognize ours as an experimentally verifi- 
able conjecture of what might happen at the lipid com- 
ponent level since the cited experimental evidence only 
corresponds to chemical processes involving much larger 
molecules (proteins). 





FIG. 1: Schematic representation of the model. Placed in 
one of the layers of the membrane, the two constituents con- 
stituents A and B induce positive and negative curvatures, 
respectively (dashed lines). The species are interconverted 
by means of a nonequilibrium reaction: k+ and fc_ are the 
forward and backward rate constants, respectively. The me- 
chanical effect of the reactive process on the membrane shape 
is illustrated by the hollow arrows. 



We write Eqs. (|2a[l and l|2bjl in detail in terms of di- 
mensionless parameters where the energy is measured in 
units of the thermal energy fcgT, the time in units of r^, 
and length in units of y 'Dt^: 



90 
dt 

dh 
~dt 



{kHI - a)V 2 </> -I- 3/30 2 V 2 + 6/?0|V0| 2 

-7V 4 </>-Ki7 o V 4 /i-r(< ? !)-0 o ), (3a) 
-kV 4 H kH V 2 (I) + T^ - <p )C (3b) 



In the next two sections we show analytically and con- 
firm numerically that the new force term can have pro- 
found effects on the instability and pattern formation 
properties of the model. In particular, this term may 
lead to dynamical patterns, whereas in its absence only 
stationary patterns are possible. 



III. ANALYTICAL TREATMENTS 

The stationary uniform state corresponds to <p — <Ao 
and arbitrary h. The linear stability of these uniform 
solutions is tested by adding small plane- wave perturba- 
tions of wave numer q to the uniform state and linearizing 
Eqs. (|3a|) and l|3b|) in these perturbations. This proce- 
dure leads to the 2x2 linearization matrix C with the 



following coefficients, 

C 21 = - K H Q q 2 +Tt 
L22 = ~nq . 



+ 3(3$)+ iq 2 ] -r 



(4) 



The eigenvalues ui q of the Jacobian associated with the 
matrix C correspond to the linear growth rates of the 
perturbations. 

The solutions of the eigenvalue problem are given 

by oj q = \ [Tr\L\ ± V^Zj) , where A[£] = Tr[C} 2 - 

4Det[C]. At the instability boundary, TLe(u> q ) vanishes 
for one finite wave number that is defined as the first un- 
stable mode. If the imaginary part of u> q is not zero at this 
wave number, we have a wave bifurcation. The condition 
for this bifurcation is obtained by requiring Tr[£] = 
and A[£] < 0. If the imaginary part of the growth rate 
is zero, one has a Turing-like bifurcation. The condition 
for this bifurcation is Det[C] = 0. 

In the absence of the interconversion reaction it is 
of course well known that immiscibility leads to com- 
plete phase separation. This occurs because a continuous 
range of modes starting from q = is unstable, and phase 
separation does not stop until there is complete segrega- 
tion into two large domains. In our previous work |2fJ 
we showed that the interconversion reaction provides a 
mixing mechanism that stabilizes the longest wavelength 
modes. If the reaction rate parameter T is sufficiently 
large, 



r > Ti 



ml) 



(5) 



then mixing is complete, the instability disappears, and 
the bilayer becomes stable and essentially flat. If T < Ti, 
the reaction no longer stabilizes all modes, but it does 
stabilize the longest wave vector modes. The first mode 
to become unstable is 



a - 3/3$ 



(6) 



The longest wavelength unstable mode now lies at the 
finite value 



3(3$ 



20 



3(3$ 



(7) 



This means that the phase separation process stops when 
the separated regimes are of characteristic size ~ <7min- 
The patterns are Turing-like (stationary) because the 
imaginary part of the growth rate is zero at the bifur- 
cation point. A typical phase diagram found earlier for 
this case is shown in Fig. [21 Note that the pattern size is 
independent of curvature parameters. Curvature reduces 
the unstable mode growth rates (see the £ = curves 
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in Fig. |3J), but without changing either the characteris- 
tic size of the patterns nor the marginal condition JSJ. 
The effect of curvature is thus restricted to the kinetics 
of the phase sep aration process, as has been extensively 
analyzed in [29j. 
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next section and exhibit two different views of a typ- 
ical phase diagram in Fig. 0] Regime 111 is the sta- 
ble regime where the mixing due to the interconversion 
is sufficiently strong to produce a homogeneous phase. 
Here "sufficiently strong" depends on the value of the 
parameter £. Beyond this stable regime, there are now 
two instability regimes. One, called region II in the di- 
agrams, is the Turing instability regime leading to sta- 
tionary nonequilibrium patterns, that is, the same sorts 
of patterns observed in the absence of the reaction-shape 
coupling (upper region of Fig. |2J) 29] . The condition for 
the Turing-like bifurcation, Det[C] — 0, obtained from 
the linear stability conditions now is T < T2, with 



r 2 = 



(8) 



This is the curve bounding region II in the figure. The 
first Turing-like unstable mode occurs at 
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(9) 



FIG. 2: Phase diagram for the case of immiscible components 
in the plane (a, V) for 7=1 and 0o = when there is no 
force on the membrane caused by the interconversion reaction. 
From El. 
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FIG. 3: Dispersion relation functions u) q for different values 
of r and £. The other parameters are held fixed at a — 7 = 1, 
cj>o — 0, k = 0.5 and Hq = 0.2. Only the case with £ = 5 has 
a nonzero imaginary part of the dispersion relation (plotted 
with symbols). 



The new questions of interest here are the consequences 
of the reaction-shape coupling in the kinetic equations. 
To place the analytic results obtained in this section 
in context, we anticipate the numerical results of the 



The new regime, region I, is one in which one observes 
dynamical (i.e., time-dependent) patterns. Provided the 
reaction parameter is sufficiently small and the reaction- 
shape coupling parameter is sufficiently large, the linear 
stability analysis in this regime leads to wave bifurcation 
solutions, Tr[C] — and A[£] < 0. The explicit condi- 
tions for these time dependent solutions obtained from 
the linear stability conditions are T < Tq, with 



r n = 



ml 



4 (7 + k) 



and £ > £oj with 



(10) 



(11) 



The first unstable mode in the above expression now is 



2(7 



7' 



(12) 



The curves bounding region I in Fig. 0] arise from these 
expressions. A typical dispersion relation function ui(q) 
in this regime is shown in Fig. [3] The function now has a 
nonzero imaginary part as well as a real part, as appropri- 
ate for spatiotemporal pattern formation. Typical spa- 
tiotemporal field configurations in the unstable regions 
of the phase diagram can only be found numerically. We 
do so in the next section. 

Before doing so, however, we can go a step further 
in the analytic approach to the problem by introduc- 
ing a weakly nonlinear analysis based on the amplitude 
equation formalism. Such a nonlinear expansion is valid 
near the bifurcation threshold to pattern formation and 
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FIG. 4: Phase diagrams in the planes (a, V) for £ = 5, and for a = 1. The other parameter values are 7 = 1, k = 0.5, 

Ho — 0.2 and 4>o = 0. The boundaries in these diagrams have been calculated numerically (see Sec. lIVI . 



sheds light on the expected pattern configurations and 
their dynamical properties. For simplicity, we restrict 
the calculation of the amplitude equations to one spatial 
dimension. Nonetheless, we can readily deduce the pos- 
sible arrangements in two dimensions by an appropriate 
interpretion of the coupling terms. 

The details of the derivation of the amplitude equa- 
tions are presented in the Appendix. The analysis re- 
veals that a solution for equations (|3a[l and (|3b|) near the 
bifurcation to pattern formation reads 

4>(x,t) = (f> + A + (x,t)eKp[i(q x + u;ot)] 

+A-(x, t) exp [i(qox — toot)] + c.c, (13a) 

h(x,t) = B + (x, t) exp [i(qox + ujQt)} 

+B-{x, t) exp [i(qox — u>ot)] + c.c, (13b) 

where c.c. stands for the complex conjugate, ujq = 
|Im(w go )|, and the amplitudes A±(x,t) and B±(x,t) sat- 
isfy the equations 

dA± 3/3g 2 (r -r) 
" I- 1 - - 1 ) A ± p 



dt 



[ s 24^ggg(8K g g + w ) ^ A± | A± |2 



Cl 



-c 2 A ± \A T \ 2 

-(3/3^ + H 2 k + 6 7 g 2 - a) 



d 2 A ± 
dx 2 



-6H Q Kq Q 



2 d 2 B ± 



dB± 
dt 



dx 2 ' 



= -£(r -r)A± 

TT d 2 A± „ 2 8 2 B ± 



(14a) 



(14b) 



where 

ci = 8^ 4 [r (l + ff + 4^(3/3^ + 47<? 2 -c0] - 2uj 2 
-two [r + Mo (3/30g + H 2 k -a + 4q 2 ( 1 + «))) , 

(15a) 

6/3g 2 (r - r) 



C-2 



x 1- 



r (l + HoO + Aq 2 {Zp4>l + 479o 2 -a) J' 

(15b) 

Note that no pattern develops if T > Tq since the ampli- 
tudes then decay with time, A±, B± — * 0. Note also the 
coupling term between the rightward, A_, and leftward, 
A+, traveling waves. This interaction is responsible for 
generating either standing or traveling waves in the mem- 
brane dynamics. Standing waves appear as a result of a 
kink- antikink dynamics where both waves have the same 
amplitude, A_ = A + . Traveling waves require these am- 
plitudes to be different. Therefore, if the coefficient c 2 of 
this interaction term is small (strictly speaking, zero) we 
expect the two amplitudes to be equal, so that 

</>(x,t) — (f>o + 4:A + (x,t) cos(q x) cos(ujot), (16a) 
h(x,t) — 4B + (x,t) cos(qox) cos(ujQt), (16b) 

that is, the solution is a standing wave. On the other 
hand, if c 2 is not zero, a traveling wave appears. Ob- 
serve that if $0 = 0, that is, if the composition of the 
membrane is 50% lipid A and 50% lipid B, then the 
A + A- interaction term is always relevant and no 
standing waves are possible. This scenario is confirmed 
by the numerical simulations described in the next sec- 
tion (cf. Figs. [S] and where we see a composition- 
curvature traveling wave in the membrane. For the pa- 
rameters in that figure, c 2 f=a 0.25). On the other hand, 
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this restriction does not apply if the membrane compo- 
sition is such that ^ / (cf. Figs. [7| and 0D where 
<f>o = 0.14 and C2 ~ 10 -2 ; the simulations show that in 
those cases there is a standing oscillatory contribution to 
the pattern). As for the different spatial configurations, 
one obtains either rolls or hexagons depending on the 
symmetries of the system. Thus, note that in Eqs. i|3a|) 
and Ij3bfl the inversion symmetry 

<j) ->-</>, h^-h (17) 

is satisfied only if <fio — and in that case roll-like ar- 
rangements are obtained. On the other hand, if </>o 
there is no inversion symmetry, and the resulting struc- 
tures are expected to be hexagonal. 

IV. NUMERICAL RESULTS 

Representative numerical results are presented in this 
section, showing good agreement with the predictions 
of the linear stability and amplitude equation analy- 
ses. We numerically solve Eqs. l|5aj) and (|3b|) in a two- 
dimensional square lattice of 100 x 100 sites with a mesh 
size Ax = 1 and periodic boundary conditions. The 
spatial derivatives are calculated by means of a sim- 
ple centered scheme, and a first-order Euler algorithm 
with time step At = 10~ 4 is used for the temporal 
integration. These coordinate and time steps insure 
good numerical accuracy. Simulations are started from 
a slightly randomly perturbed homogeneous distribution 
4>(r) = cj) = cf)Q and h(r) = h = 0. In this section we 
present numerical results that correspond to highly im- 
miscible components (deep quench, a = 1 and [3=1, 
leading to an equilibrium value of 4> eq = ±1), and an in- 
terface thickness of the order of the space discretization 
(C = 1, which leads to 7 = 1). The bending rigidity 
modulus is taken equal to 0.5 (in units of fcgT), and very 
different shapes are implemented for the two constituents 
by setting Hq = 0.2. 

When dealing with mesosopic or coarse-grained mod- 
eling schemes one has to be concerned about the applica- 
bility of the results. The limit on the predictive power of 
these models is mainly determined by the identification 
of the correct time, length and energy scales accessible to 
the experiments. In our case, we have adimensionalized 
the kinetic equations in order to have D = A = fc^T = 1, 
so we can return dimensions to these equations by choos- 
ing typical values for these parameters. The diffusion co- 
eficient is known to be in the range 10 -7 — 10 _8 cm 2 /s (for 
lipids in a liquid-disordered phase membrane ^i|). In the 
absence of the nonlinear and mechanical effects consid- 
ered in this paper, solvent hydrodynamic effects (Zimm 
dynamics) |3J can be incorporated with a renormalizcd 
height mobility A = (4/igr) in Eq. 12bll in Fourier 
space |32] . Here /1 stands for the solvent kinematic viscos- 
ity, which is lep for water at 20°C. For n = 0.5/cbT the 
linear stability analysis shows typical unstable modes at 
q ss q ss 0.5, which corresponds to a pattern size in the 



range of 2 — 20/im. Such a size may be accessible in ex- 
periments on giant vesicles and planar bilayers. However, 
when mechanical effects are included, a more elaborate 
analysis such as that provided by dynamical renormal- 
ization methods [32j is in order. A simple and plausible 
hypothesis might be to include such effects by simply 
renormalizing both A and £ in Eq. (|2b|) , both with a 1/q 
decay in Fourier space since one might expect a non-local 
viscous kernel to decay as a power law with distance [33| . 
In any case, a more detailed consideration of this issue is 
beyond the scope of this work. 

However, in order to accurately delimit the range of 
applicability of our model, it is useful to add some com- 
ments regarding the value of n. The elastic proper- 
ties, and thus the bending rigidity, of lipid bilayers are 
strongly determined by the size, shape and molecular 
elasticity of their constituents. Membranes composed of 
phospholipids (which have two hydrophobic chains) are 
characterized by a bending rigidity of the order of tens of 
ksT [35| . Cell membranes containing large molar frac- 
tions of cholesterol (rather rigid molecules) display even 
higher rigidities 16]. In our model, the typical sizes of the 
predicted dynamical patterns emerging from a wave in- 
stability for such large values of k exceed the experimen- 
tally accessible sizes. Moreover, for k > 10, wave-like un- 
stable modes are observed to be supressed at long times 
by the other existing Turing-like unstable modes, so that 
only stationary structures are asymptotically generated 
in this parameter region. However, single- and/or short- 
chain lipid surfactants are known to form much more 
flexible bilayers, with bending rigidities of the order of 
fcfiT or even less [3(|. Furthermore, recent theoretical 
models [37| show how bilayers composed of surfactants 
of rather distinct shapes (as in our case) may exhibit 
smaller rigidities than one-component membranes. This, 
therefore, delimits the range of applicability of the model 
results in this paper. 

Our quantitative assigmcnt of £ has been made with 
the energetic feasibility of a number of external sources 
in mind. To see this, note that the characteristic time 
for a single reactive event is of order 1/r. During this 
time interval the change in the height of the membrane 
due to the contribution of the reactive term is of order 
£ [see. Eq. (|2EI); (f> - <t>o is of O(l)]. Moreover, the lo- 
cal energetic "cost" of such a height increase due to the 
reactive term is of order 2k£ 2 / (Ax) 2 [see Eq. (QJ]. In 
most of our simulations we have used £ = 5 (in simu- 
lation length units), which, for k = 0.5, corresponds to 
an energy cost of order 25 (in units of ksT). The de- 
phosphorylation of an ATP releases about 50kJ/mol or 
10 _19 J/ATP, which at room temperature is precisely of 
order 25fcsT. If the energy source is light, the power re- 
quired to produce 25k bT within a time interval of order 
r _1 for r = 0.14 as used in our simulations is only ap- 
proximately 0.3 x 10~ 6 mW/cm 2 , which is much lower 
than the power produced by typical commercial light 
sources used, for example, in photosensitive Langmuir 
monolayer experiments. While some of the energy pro- 





FIG. 5: Spatiotemporal behavior for the case of an unstable membrane in region I of the pase diagram Fig.0]and the following 
set of parameters: a = /3 = 7 = 1, <j>o = 0, k = 0.5, Ho = 0.2, F = 0.14 and f = 5 (dashed curve in Fig. 0. The two 
upper panels correspond to the <j>- an d /i-field distributions at t = 15000, once the structures are robust. The stripes travel 
from left to right at constant velocity. In the third panel we plot the temporal evolution (horizontal axis) of a one-dimensional 
cross-section ^-profile (vertical axis) at y — 50, starting at t — 15000 until t — 19000. In the fourth panel, the corresponding 
ft-profile (horizontal axis) is plotted against time (downward below sheet axis), starting at t = 15000 until t — 15500. Darker 
(lighter) regions are richer in the A (B) lipid in the concentration snapshots, and some exaggeration along the vertical direction 
has been applied to the height plots. 



vided by these external sources (ATP, light, etc.) would 
be used for processes other than the purely elastic motion 
of the membrane, including dissipation into the thermal 
surroundings and the movement of other masses, this or- 
der of magnitude estimate shows that the mechanism is 
feasible and provides a basis for the value of the param- 
eter £ chosen for our simulations. 

For a highly flexible bilayer, typical phase diagrams 
are those of Fig. 0] As noted earlier, region I corre- 
sponds to a regime where at least one of the unstable 
modes has nonzero imaginary part, so both types of in- 
stablity are present. Region II is the unstable phase also 
present in Fig. [21 Here, the unstable modes have no 
imaginary parts, so that stationary finite-sized patterns 
are predicted. In region III, there are no unstable modes 
(stable phase). Here we have exhibit spatiotemporal pat- 
terns associated with the new region I. 

The nature of the spatiotemporal patterns is deter- 
mined by the value of (f>o and by the associated relative 
magnitudes of the amplitudes of waves traveling in differ- 
ent directions. These assertions are supported by the am- 
plitude equation analysis of the previous section. Thus, 
for a critical quench (</>o = 0), the system typically dis- 
plays some transients with domains that travel in differ- 



ent directions until they organize into a coherent train of 
traveling stripes. The post-transient spatiotemporal be- 
havior is described in the different panels in Fig. [S] We 
plot the concentration and height fields (upper panels), 
and one-dimensional cross-sections showing the tempo- 
ral evolution, for both order parameters. The traveling 
waves are consistent with the predictions of the ampli- 
tude equation analysis. Numerical profiles at a given 
post-transient stage reveal the mechanism that leads to 
the motion of the generated spatial structures. In Fig. 
we observe that <f>- and ^.-profiles are slightly displaced, 
the field <f> being ahead in the direction of propagation. 
That there is such a phase difference between the two 
fields is consistent with the appearance of the contribu- 
tion with an imaginary coefficient on the right hand side 
of Eq. JHEJl. 

Off-critical quenches (</>o = —0.14) display different 
spatial and temporal behaviors, again consistent with the 
amplitude equation results. As before, we first observe 
transients, but now involving concentration droplets and 
bud-like surface deformations. The fields settle into an 
oscillating pattern of buds rich in the minority species 
and whose geometry depends on the value of the param- 
eter £ that measures the mechanical effect of the reaction 
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model. Azobenzene derivatives, which are also known to 
display dynamical organization in noncquilibrium Lang- 
muir monolayers [38| , may be suitable compounds to test 
our predictions since their shapes can be strongly modi- 
fied by means of well-known (externally controlled) pho- 
toisomerization reactions [3j|, and the resulting isomers 
normally exhibit phase separation. 

Finally, we note that it would be interesting to study 
the crossover between Rouse and Zimm dynamics in our 
model system so as to elucidate how the different pa- 
rameters are renormalized due to solvent hydrodynamic 
effects, and to clarify the resulting consequences for the 
pattern formation mechanism proposed herein. Such 
a study has been carried out for tethered membranes, 
where different coarsening exponents are found [32, H3 • 



FIG. 6: <j>- and /(.-profiles for the case in Fig.|^]at a given time 
t = 18000. Fhe arrow represents the direction of propagation 
of the traveling pattern. 



on the shape of the membrane. This behavior is shown 
in Fig. For £ — 5 a spatial Fourier transform of the 
pattern indicates clear square symmetry. For £ = 3, the 
domains are rather hexagonal and they oscillate as well 
as move in one direction, as shown in Fig. [SJ 
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CONCLUSIONS 



APPENDIX A: DERIVATION OF THE 
AMPLITUDE EQUATIONS 



In a previous model for a flexible reactive bilayer com- 
posed of two differently shaped lipids, we demonstrated 
the occurrence of stationary finite-sized domains of com- 
position and curvature as a result of the competition be- 
tween thermodynamic phase segregation and a nonequi- 
librium reaction [2^| . A generalization of the model that 
includes the local effect of the reaction on the membrane 
shape has been presented here. We have shown how the 
mechanical influence of the reactive process on the mem- 
brane may lead to the formation of spatiotemporal struc- 
tures. The linear stability analysis of the model equations 
shows the existence of a wave instability for sufficiently 
large reaction-shape coupling. A weakly nonlinear anal- 
ysis based on the amplitude equation formalism provides 
insight into the spatial and temporal symmetries of the 
emerging patterns. Correspondingly, numerical simu- 
lations display the spontaneous generation of traveling 
lamelar phases, and oscillating or moving buds, showing 
an important potential for the formation of spatiotempo- 
ral patterns in deformable nonequilibrium reactive lipid 
membranes. 

This extension of the model should be considered as a 
formal issue in reactive deformable surfaces. However, we 
believe that experimental work on giant vesicles and on 
planar membranes can be designed to qualitatively cap- 
ture the spatiotemporal phenomenology obtained in our 



For convenience of notation, we rewrite the one dimen- 
sional version of the model equations l|3al) and l|3bjl for 
reactive membranes in terms of the shifted held variable 
if — (f> — (fio and the control parameter e = (ro — T) /Tq 
that accounts for the "distance" to the threshold of the 
bifurcation to pattern formation so that a pattern devel- 
ops if e > 0: 



(kHq - a)ip xx + 3f3(ip + (f>o)ip x 
+6/% + O ) (<p 
-nH h xxxx - r (l - s)(f, 



(Ala) 



Kh xxxx + kHq 

Vxx + £r (l - e)<P- (Alb) 



The subscripts indicate partial derivatives. The homoge- 
neous state is (p = and arbitrary h — h. With no loss 
of generality we take h = 0. 

Slightly above the bifurcation to p attern formation, 
the following expansion holds (see |29l |41| and references 
therein): 



n=l 

CO 

h = J2^ /2 h {n) , 



(A2a) 
(A2b) 
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FIG. 7: Spatiotemporal behavior for the same case as in Fig. 03 except for 0o = —0.14 (still in region I). The snapshots 
correspond to t — 22500 and the temporal profile evolutions go from this time up to t — 26500 for <j>, and to t = 23000 for h. 
The same panel organization and specifications as in Fig. apply. Bud-like curvature domains oscillate. 



and a separation of spatial scales can be implemented 
between the most unstable mode (fast growth), qo, and 
the rest of the modes within the unstable band (slower 
growth). In terms of the control parameter, we let 
X = £ x / 2 x denote the spatial modulation scale of the 
slow modes and T = et the associated time scale. The 
separation of scales between the fastest growing mode 
and the slower modes can be implemented in Eqs. (|Ala|) 
and l)Alb|) by expanding the spatial and temporal deriva- 
tives according to the chain rule so that d x — * d x +e 1 ' 2 dx 
and dt — > dt + sdx- Implementation of this scale separa- 
tion and substitution of Eqs. l|A2a(l and (| A2b|) into IjAlal) 
(|Alb[) leads to a rather cumbersome expansion in terms 
of e. The lowest order contribution is of order e 1 / 2 , and 
balancing all terms of this order gives 



inition of the linear operator L with elements 



(3/?0g + H 2 k- a) 



dx 2 



d 4 



L 12 - -HoKq^ 
L 2 i = £r + H k 



->22 — —K 



dx 4 



dx 2 
d_ 
dt' 



d_ 

dt 



Eqs. i|A3a(l and (|A3b|l can trivially be written as L%i = 0, 
where (x„) T = (</? (n) , /i (n) ). 

The next order of the expansion gives us terms of 
0{e) and yields 1L X 2 = V>2 where = 

(4 a) ,4 6) ),and 



(3j9^+fZ?K-a)^-7^ 
- H Kh£ a 



(i) 



XX 



^ = o, 

xxxx 



W = o. 



(A3a) 
(A3b) 



4 a) 



4 b) = 2K[-H <p% + 2h% xX 



-6/300 



+2 (a - ml ~ H 2 k) <pW + 4 7 ^ 



(1) 



X 



Equations (|A3a|) and (|A3b|) are in fact simply the lin- 
earized versions of Eqs. i|ATaJl and (|A"Tb|) . With the def- 



At order e 3/2 we get L% 3 = V>3 ({<^ (1) , t^ (2) ; /i (2) }), 
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FIG. 8: Spatiotemporal behavior for the same case in Fig. 03 except for 0o = —0.14 and £ = 3 (still in region I). In the 
upper panels, the droplet-like structures move from the upper-left corner to the lower-right one at constant velocity, while they 
oscillate. The snapshots correspond to t = 18000 and the temporal profile evolutions go from this time up to t = 22000 for <f>, 
and to t = 18500 for h. The same panel organization and specifications as in Fig. |S] apply. 



where the components of ip^ 



4 a) 



4 b) 



^3°\'03 >) ) are given by equation, such that at order e n/2 

L X „ = rp n . . . , ^(-D; hV, . . . , h^-V}) . (A4) 



3/3 (^+2^) 
27 (3^2xx 



2^ (2) x 

t xxxX 



+-2#, 



txXX 



(1) 



2h 



-H, 



(2) 

xxxX 



xx 



2h 



(2) 



cX 



2^1 



We could continue up to any order with the expan- 
sion. However, at order e 3 / 2 we are able to extract a 
closed evolution equation for the amplitudes of the pat- 
tern, as shown below, and we therefore stop at this order. 
Note that beyond order e 1 / 2 , which corresponds to the 
linear problem, in all cases we obtain a non-homogeneous 



At order e 1 / 2 the problem is homogeneous and, with ap- 
propriate boundary conditions, we can write the solution 

as 



(1) 



A 1+ {X,T)e li - qaX+u>ot) 

+A 1 -{X 1 Ty {qoX ~ UJot "> 

B 1+ (X,T)e i( - qox+Uoi ^ 



c.c, 



-B x _{X,T)e 



i(q a x-ui t) 



+ C.C. 



(A5a) 
(A5b) 



where c.c. stands for the complex conjugate and u>o = 
|Im(u; go )|. However, the amplitudes A\± and Bi± are 
undetermined at this point. As mentioned above, the 
subsequent orders are no longer homogeneous and there- 
fore the existence of an analytical solution can not be en- 
sured. However, we can enforce solvability by evoking the 
Fredholm alternative theorem |40j . In our case the appli- 
cation of the theorem simply states, as a recipe, that for 
Eas. l)A4|l to have a solution, the functions tp n can not con- 
tain the fundamental mode exp [± (iqox ± u>ot)]. In other 
words, the particular solutions of the non- homogeneous 
problem are orthogonal to the solutions of the homoge- 
neous problem. Thus, by substituting the solutions (|A5a(l 
and (|A5b|) into the next order of the hierarchy and im- 
posing the solvability condition we obtain the following 
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particular solution, 

tpW = c llp [A 1+ {X,T)f e 2 «(?o^Dt) 
+c 2lfi [Ai- (X, T)] 2 e 2 K<iox-uot) 
+c 3v A 1+ (X,T)A 1 _(X,T)e 2l ' l0X 

h {2) = cm \{A 1+ {X,T)) 2 e 2l(qoX+Uat) 

+ (A 1 _(X,T)) 2 e w(»»-'-b*) 
+ C2 ^ 1+ (X,T)A 1 _(X,r)e 2 ^ a; 



+ c.c. 



where 

num 
den 

Clh 

C2h 



+ C.C., 



r .den 

[U(3(poqt (8nq* + iwo)] 
[8H oK q%(4HoKq 2 -T Q t) 
-{8nq* + iluo) (T + 4 g 2 (3/3 ( /> 2 + H 2 k 
+4 7g 2 -a) + 2iu )] , 
[l2f3<f> Q q 2 (8Kq* - imp)] 
c 



r„(l + HoO + 4g 2 (3/3^ + A iq l - a) ' 
[6^ ogo 2 (-^r + 4g oK g 2 )] 

c* 

3^ ogo 2 (gr - AHpngp 

2 K ql [r (l + H Q Z) + 4g 2 (3/30 2 + 4 7g 2 - a)] ' 
8 K q 4 Q [r (l + #oO + 4g o 2 (3/30 2 + 47gg - a)] 
-2wq - icj [To 

+4g 2 (3/3^ 2 + i/ 2 « - a + 4 9 2 ( 7 + «))], 



and c* is the complex conjugate of c. The values of A\± 
and Si± can not be determined at this order either. How- 
ever, at order e 3 / 2 the Fredholm theorem provides closed 
equations for the conditions that must be satisfied by 
these amplitudes. These conditions are the amplitude 
equations for our pattern forming system, and are given 
in terms of the original variables x and t in Eqs. l|13a|) 
and (|13b(l . In those equations the "1" in the subscript 
has been dropped for economy of notation. 
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